%-------------------------------------------------------------------------;
% Horizon Bias and the Term Structure of Equity Returns;
% This code replicates Table 8;
%-------------------------------------------------------------------------;

clear all;
clc;
close all;

%-----------------------------;
% Import term structure data;
%-----------------------------;

filename='Term_structure_variables';
data_ts=xlsread(filename);

b=size(data_ts,1); %Length of the data;

Div_ret=data_ts(2:end,6); %BBK extended;
SP_ret=data_ts(2:end,4);

Div_ret=log(1+Div_ret);
SP_ret=log(1+SP_ret);
                                                 
etp=SP_ret-Div_ret;

SP_price=data_ts(1:end,1);
SP_Div=data_ts(1:end,3);
SP_strip_18=data_ts(1:end,5);

SP_pd=log(SP_price./SP_Div);
SP_pd_strip=log(SP_strip_18./SP_Div);

SP_pd=SP_pd(1:b-12);
SP_pd_strip=SP_pd_strip(1:b-12);

%-------------------------------;
% Import explanatory variables;
%-------------------------------;

filename='Explanatory_variables';
data_ev=xlsread(filename);

data_ev=data_ev(1:b-12,:);

year=data_ev(1:end,1);
month=data_ev(1:end,2);

Ltg=data_ev(1:end,15);   
Stg=data_ev(1:end,18);   

%Take logs;
Ltg=log(1+Ltg);
Stg=log(1+Stg);

%Rational forecasts;
Ltg_r=data_ev(1:end,16);
Stg_r=data_ev(1:end,19);

Ltg_r=log(1+Ltg_r);
Stg_r=log(1+Stg_r);

%Bias in forecasts;
Ltg_b=data_ev(1:end,17);
Stg_b=data_ev(1:end,20);

Ltg_b=log(1+Ltg_b);
Stg_b=log(1+Stg_b);

%Skewness;
skew1=data_ev(1:end,7);
skew1=(skew1-mean(skew1))./std(skew1);

skew5=data_ev(1:end,8);
skew5=(skew5-mean(skew5))./std(skew5);
Diff_skew=skew5-skew1;

%Uncertainty;
uncert1=data_ev(1:end,9);
uncert1=(uncert1-mean(uncert1))./std(uncert1);

uncert5=data_ev(1:end,10);
uncert5=(uncert5-mean(uncert5))./std(uncert5);
Diff_uncert=uncert5-uncert1;

%Other variables capturing business cycle;
rec=data_ev(1:end,11);
default=data_ev(1:end,12);
term=data_ev(1:end,13);
cay=data_ev(1:end,14);

%Horizon bias;
acorr=autocorr(Stg);
theta_0=mean(Stg);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg;

Stg_impl=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB=Ltg-Stg_impl;

%Horizon bias - rational;
acorr=autocorr(Stg_r);
theta_0=mean(Stg_r);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_r;

Stg_impl_r=(Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two);

HB_r=Ltg_r-Stg_impl_r;

%Horizon bias - bias;
acorr=autocorr(Stg_b);
theta_0=mean(Stg_b);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_b;

Stg_impl_b=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB_b=Ltg_b-Stg_impl_b;

%-------------------------;
% Predictive regressions;
%-------------------------;

h=12;  

%Number of lags for overlapping observations;
lags=h-1;

%Sum returns over h periods;
ret_SP=filter(ones(h,1),1,SP_ret);
ret_Div=filter(ones(h,1),1,Div_ret);

ret_SP=ret_SP(h:end);
ret_Div=ret_Div(h:end);

ret_SP=ret_SP(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
ret_Div=ret_Div(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

SP_pd=SP_pd(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
SP_pd_strip=SP_pd_strip(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

Ltg_h=Ltg(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
Stg_h=Stg(1:min(size(Stg,1),(size(Stg,1)-h+12)));

HB=HB(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
HB_r=HB_r(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
HB_b=HB_b(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

rec=rec(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
default=default(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
term=term(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
cay=cay(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

skew1=skew1(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
skew5=skew5(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

%Dependent variable;
y=(ret_SP-ret_Div);

%------------------------------------------------;
%Table 8;
%------------------------------------------------;

%Column 1: HB rational + controls + error + NBER;

x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB_b SP_pd term default cay error rec];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table8(1:14,1)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1);b(8,1);tb(8,1)];
Table8(15:16,1)=[length(error);AR2];

%Column 2: HB bias + controls + error + NBER;

x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB_r SP_pd term default cay error rec];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table8(1:14,2)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1);b(8,1);tb(8,1)];
Table8(15:16,2)=[length(error);AR2]
